jupyter notebook companion to interactive visualization
If you are reading this, it is very likely that you are a student. Actually, Imperial College Student, and you are undertaking Linear Algebra course this year or want to refresh fading memories.
This notebook is intended to accompany an interactive visualization for the Linear Algebra course, with more detailed descriptions of what is what and how the visualization works in the first place.
You can compare this Jupyter notebook to Half-Blood Prince's notes from Harry Potter fantasy series. The code you will find here might not exactly copy what your lecturer has prescribed you in the notes, but we hope this might give you deeper insight into your study. Here we will not only attempt to make sense of the notes, but also highlight certain problematic concepts in practice, in Python 3. We really want you to have intuition of what's going on.
The most important import, typically, is NumPy. This will allow us to perform many operations that are tedious to code by oneself, such as cross product, dot product etc.
Second import is a package written for this visualization for you to use. It is called pointslinesplanes and implements classes of Points, Lines and Planes. Do not bother about it too much; all its use will be appropriately commented. Just remember that the shorthand for the name is plp
In [2]:
import pointslinesplanes as plp
import numpy as np
Because we want to see what we do, we need some plotting package. Throughout this notebook, we use plotly to generate the relevant graphs. Feel free to acquaint yourself with it, as it might be helpful in the future. However, this is not necessary as all the code that uses plotly will have appropariate comments in this notebook.
In [1]:
import plotly.offline as py # imports plotly for the offline use
import plotly.graph_objs as go # imports plotly's graphical objects (such as scatters, surface graphs etc)
py.init_notebook_mode(connected=True) # makes plotly plots interactive
In [3]:
def intersect(obj1, obj2):
"""Intersection between two objects that can be Point, Line or Plane.
Returns Point, Line or Plane or throws ArithmeticError when no intersection can be found.
@author Nick Metelski
@since 26.07.17"""
#plane-plane
if isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Plane):
return _plane_plane_intersect(obj1,obj2)
#plane-line
elif isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Line):
return _plane_line_intersect(obj1,obj2)
elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Plane):
return _plane_line_intersect(obj2,obj1)
#plane-point
elif isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Point):
return _plane_point_intersect(obj1,obj2)
elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Plane):
return _plane_point_intersect(obj2,obj1)
#line-line
elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Line):
return _line_line_intersect(obj1,obj2)
#line-point
elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Point):
return _line_point_intersect(obj1,obj2)
elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Line):
return _line_point_intersect(obj2,obj1)
#point-point
elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Point):
return _point_point_intersect(obj1,obj2)
#wrong params
else:
raise TypeError("Invalid parameter types - please pass intersections.Point, intersections.Line or intersections.Plane.")
def _plane_plane_intersect(p1,p2):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
plane-plane intersection submethod. Raises exceptions or returns Line.
@author Nick Metelski
@since 26.07.17"""
#TODO make it return plane is planes are identical
#cross-product
cross = np.cross(p1.normal,p2.normal)
if np.all(cross==0):
#planes are parallel or overlap
#if they overlap then the offset of one plane lies on the other
try:
intersection = intersect(Point(p1.offset), p2)
return Plane(p1.normal,p1.offset) #return a new plane same a the overlapping ones
except ArithmeticError:
raise ArithmeticError("Planes are parallel, do not overlap.")
else:
#sample point: x=0
#NOTE: linalg.solve can raise np.linalg.LinAlgError exception when x=0 results in singular matrix
#we have to check some other cases
#premise: the line has to intersect at least one of the planes: xy, yz or xz.
mat = [[p1.normal[1],p1.normal[2]],[p2.normal[1], p2.normal[2]]]
axis = 0 #keep track of which axis we zero out; 0=x, 1=y, 2=z
if np.linalg.matrix_rank(mat) == 1:
mat = [[p1.normal[0],p1.normal[2]],[p2.normal[0], p2.normal[2]]]
axis = 1
if np.linalg.matrix_rank(mat) == 1:
mat = [[p1.normal[0],p1.normal[1]],[p1.normal[0], p2.normal[1]]]
axis = 2
rhs = [np.dot(p1.normal,p1.offset),np.dot(p2.normal,p2.offset)]
sol = np.linalg.solve(mat,rhs)
if axis == 0:
return plp.Line(cross, [0,sol[0],sol[1]])
if axis == 1:
return plp.Line(cross, [sol[0],0,sol[1]])
if axis == 2:
return plp.Line(cross, [sol[0],sol[1],0])
def _plane_line_intersect(plane,line):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
Intersection of a line an a plane. Raises exceptions or returns Point.
@author Nick Metelski
@since 26.07.17"""
check = np.dot(plane.normal,line.vec)
if np.allclose(check,[0.0]):
#plane and line are parallel or overlap
raise ArithmeticError("Plane and line are parallel.")
else:
#there is an explicit formula for parameter of the line for which intersection is met:
#$$t = \frac{\vec{n} \cdot (\vec{d_p}-\vec{d_v})}{\vec{n} \cdot \vec{v}}$$,
#where n is normal to plane, v is line's direction, rp is plane offset, rv is vector offset
t = np.dot(plane.normal,(plane.offset-line.offset))/np.dot(plane.normal,line.vec)
return plp.Point(line.offset+t*line.vec)
def _plane_point_intersect(plane,point):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
Intersection of a plane and point. Raises exceptions or returns Point.
@author Nick Metelski
@since 27.07.17"""
check = np.dot(plane.normal, point.pos - plane.offset)
if check == 0:
return plp.Point(point.pos)
else:
raise ArithmeticError("The point does not line on the plane.")
def _line_line_intersect(line1,line2):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
Intersection of two lines. Raises exceptions or returns Point or Line.
@author Nick Metelski
@since 27.07.17"""
#raise NotImplementedError("Line-line intersections not implemented.")
cross = np.cross(line1.vec,line2.vec) #is normalized as vectors are normalized
if np.allclose(cross,[0.0,0.0,0.0]):
#lines parallel. are they identical?
try:
intersect(line2, plp.Point(line1.offset)) #this will raise an exception if they are not identical
return plp.Line(line1.vec,line1.offset)
except ArithmeticError:
raise ArithmeticError("Lines are parallel, and are not identical.")
#lines not parallel
sample = line2.offset - line1.offset #sample vector between two points on the line
dist = np.abs(np.dot(sample,cross)) #distance between the lines
if np.allclose([dist],[0.0]):
v1 = line1.vec
v2 = line2.vec
d1 = line1.offset
d2 = line2.offset
n = plp.utils.normalize(np.cross(v2,v1))
s = d2 - d1
l = np.linalg.norm(np.cross(s,v1))
l_vec = l * np.cross(v1,n)
param = l**2/np.dot(l_vec,v2)
return plp.Point(line2.offset+param*line2.vec)
else:
raise ArithmeticError("Lines are skew, no intersection. Distance: " + str(dist))
def _line_point_intersect(line,point):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
Intersection of a line and point. Intersection is a Point.
Raises exceptions or returns Point.
@author Nick Metelski
@since 27.07.17"""
#premise: difference in offsets must be parallel to direction vector
offset_diff = point.pos - line.offset
cross = np.cross(line.vec,offset_diff) # this is for checking if there exists an intersection
if np.allclose(cross,[0.0,0.0,0.0]):
return plp.Point(point.pos) #return a point coinciding with the original point passed as param
else:
raise ArithmeticError("Point does not lie on the line.")
def _point_point_intersect(point1, point2):
"""Private function; do not use. Use intersect(obj1,obj2) instead.
Intersection of a two points. Intersection is the point if the point are the same position.
Raises exceptions or returns Point.
@author Nick Metelski
@since 27.07.17"""
if np.allclose(point1.pos,point2.pos):
#just return one of the points when they coincide
return plp.Point(point1.pos)
else:
raise ArithmeticError("Points do not coincide.")
In [4]:
plane1 = plp.Plane([1,1,1],[1,0,0])
plane2 = plp.Plane([-1.0,1,1],[0.5,0,0.5])
plane3 = plp.Plane([0.2,0.1,0.3],[0.2,-0.3,0.5])
line1 = intersect(plane1,plane2)
line2 = intersect(plane2,plane3)
line3 = intersect(plane1,plane3)
point1 = intersect(line1,line2)
In [5]:
layout = dict(
width=1000,height=1000,
showlegend=False,
font = dict(family="Verdana"),
scene = dict(
xaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
yaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
zaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
aspectmode = 'cube',
camera = dict(center=dict(x=0,y=0,z=0),eye=dict(x=0,y=0,z=3))
),
plot_bgcolor='rgb(255, 255, 255)'
)
In [6]:
objs = [plane1, plane2, plane3, line1, line2, line3, point1]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)
In [29]:
#Line(dir=Vector(0.447,0.894,0.000),off=Vector(0.120,-0.060,0.000)) and Line(dir=Vector(0.316,0.949,0.000),off=Vector(0.150,-0.050,0.000))
line2 = plp.Line([0.4472135954999579,0.8944271909999159,0],[0.120,-0.06,0])
line1 = plp.Line([0.31622776601683794,0.9486832980505138,0],[0.15,-0.05000000000000002,0])
v1 = line1.vec
v2 = line2.vec
print("v2",v2)
d1 = line1.offset
d2 = line2.offset
n = plp.utils.normalize(np.cross(v2,v1))
print(n)
s = d2 - d1
print(s)
l = np.linalg.norm(np.cross(s,v1))
print(l)
l_vec = l * np.cross(v1,n)
print("v1 x n", np.cross(v1,n))
print("l_vec ", l_vec)
print("lxl", l*l)
print("dot",np.dot(l_vec,v2))
param = (l*l)/np.dot(l_vec,v2)
print(param)
point = plp.Point(line2.offset+param*line2.vec)
In [28]:
objs = [line1,line2,point]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)
print(point.goify())
In [34]:
plane = plp.Plane([0.1,0.3,-0.2],[0.2,0.1,0.4])
line = plp.Line([0.1,0.3,0],[0.2,0.1,0])
pt = intersect(plane,line)
In [35]:
objs = [plane,line,pt]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)
print(point.goify())
In [ ]: